notebook.community

Edit and run 





In [ ]:



    


import sys
globber_path = '/Users/adrian/projects/globber/'
if globber_path not in sys.path:
    sys.path.append(globber_path)

# Third-party
from astropy.io import ascii
import astropy.coordinates as coord
import astropy.units as u
import h5py
import matplotlib as mpl
import matplotlib.pyplot as pl
import numpy as np
pl.style.use('apw-notebook')
%matplotlib inline
from scipy import interpolate
from scipy.misc import logsumexp

from astroML.utils import log_multivariate_gaussian

from globber.ngc5897 import r_t, r_c, cluster_c
from globber.ngc5897 import fiducial_DM as DM



    











In [ ]:



    


# XCov_filename = "/Users/adrian/projects/globber/data/ngc5897/XCov.h5"



    











In [ ]:



    


# with h5py.File(XCov_filename, "r") as f:
#     iso_X = f['isochrone']['X'][:]
#     clu_X = f['cluster']['X'][:]

First, compare isochrones and data





In [ ]:



    


dart_iso_filename = "/Users/adrian/projects/globber/data/ngc5897/dartmouth_iso_ps1.dat"
dart_iso = ascii.read(dart_iso_filename, header_start=8)

# pars_iso_filename = "/Users/adrian/projects/globber/data/ngc5897/ngc5897_iso_ps1.dat"
# pars_iso = ascii.read(pars_iso_filename, header_start=13)
# pars_iso[114:] = pars_iso[114:][::-1]

pars_iso_filename = "/Users/adrian/projects/globber/data/ngc5897/parsec_iso_ps1.dat"
pars_iso = ascii.read(pars_iso_filename, header_start=13)[:114]



    











In [ ]:



    


ps1 = np.load("/Users/adrian/projects/globber/data/ngc5897/PS1_stars_pv3_dered_sm.npy")
ps1_c = coord.SkyCoord(ra=ps1['ra']*u.degree, dec=ps1['dec']*u.degree)
cluster = ps1[ps1_c.separation(cluster_c) < 3*r_c]



    











In [ ]:



    


fig,axes = pl.subplots(1,2,figsize=(8,8),sharex=True,sharey=True)

axes[0].plot(cluster['dered_g']-cluster['dered_i'], cluster['dered_i'], ls='none', marker='.', alpha=0.5)
axes[0].plot(dart_iso['gP1']-dart_iso['iP1'], dart_iso['iP1']+DM, marker=None, lw=2.)
axes[0].plot(pars_iso['gP1']-pars_iso['iP1'], pars_iso['iP1']+DM, marker=None, lw=2.)

# axes[1].plot(clu_X[:,2], clu_X[:,0], ls='none', alpha=0.5)
axes[1].plot(dart_iso['gP1']-dart_iso['iP1'], dart_iso['iP1']+DM, marker=None, lw=2.)
axes[1].plot(pars_iso['gP1']-pars_iso['iP1'], pars_iso['iP1']+DM, marker=None, lw=2.)

pl.xlim(0,0.7)
pl.ylim(21,17)



    











In [ ]:



    


itrp_i = np.linspace(17., 21.5, 128)
itrp_gi = interp(itrp_i)



    











In [ ]:



    


pl.figure(figsize=(4,8))

pl.scatter(itrp_gi, itrp_i, alpha=0.5)

pl.xlim(0,1.)
pl.ylim(21,17.5)



    











In [ ]:



    


x = np.arange(0,1,0.02)
y = np.arange(17.5,21.,0.04)
shp = (y.size, x.size)
xygrid = np.vstack(list(map(np.ravel,np.meshgrid(x, y)))).T



    











In [ ]:



    


xymodel = np.vstack((itrp_gi,itrp_i)).T
xymodel.shape



    











In [ ]:



    


h = 0.02

V = np.diag([h]*2)**2
W = np.array([[1, -1],   # g-i
              [0, 1]])  # i magnitude

# each covariance C = WCW^T
# V = np.einsum('mj,jk->mk', W, V)
# V = np.einsum('lk,mk->ml', W, V)
# V



    











In [ ]:



    


ll = log_multivariate_gaussian(xygrid[:,np.newaxis], xymodel[np.newaxis], V)
ll = logsumexp(ll, axis=1)
ll.shape



    











In [ ]:



    


pl.figure(figsize=(4,8))

pl.pcolormesh(xygrid[:,0].reshape(shp), xygrid[:,1].reshape(shp), 
              np.exp(ll).reshape(shp), cmap='Blues')

pl.xlim(0,1.)
pl.ylim(21,17.5)

Check that contrast is highest near turnoff





In [ ]:



    


with h5py.File(XCov_filename, "r") as f:
    bg_X = f['control']['X'][:]
    bg_Cov = f['control']['Cov'][:]
    
    idx = (bg_X[:,0] >= 17.) & (bg_X[:,0] <= 21.5) & (bg_X[:,2] >= -0.1) & (bg_X[:,2] <= 1.1)
    bg_X = bg_X[idx]
    bg_Cov = bg_Cov[idx]
    
    bg_X = bg_X[::50,[2,0]]
    bg_Cov = bg_Cov[::50,[2,0]][:,:,[2,0]]



    











In [ ]:



    


pl.figure(figsize=(4,8))

pl.plot(bg_X[:,0], bg_X[:,1], alpha=0.4, marker='.', ls='none')

pl.xlim(0,1.)
pl.ylim(21,17.5)



    











In [ ]:



    


bg_h = 0.08
bg_V = np.diag([bg_h]*2)**2
W = np.array([[1, -1],   # g-i
              [0, 1]])  # i magnitude

# each covariance C = WCW^T
bg_V = np.einsum('mj,jk->mk', W, bg_V)
bg_V = np.einsum('lk,mk->ml', W, bg_V)
bg_V



    











In [ ]:



    


bg_X[np.newaxis].shape, xygrid[:,np.newaxis].shape, _V[np.newaxis].shape



    











In [ ]:



    


_V = bg_Cov + bg_V[np.newaxis]
bg_ll = log_multivariate_gaussian(bg_X[np.newaxis], xygrid[:,np.newaxis], _V[np.newaxis])
bg_ll = logsumexp(bg_ll, axis=1)
bg_ll.shape



    











In [ ]:



    


pl.figure(figsize=(4,8))

pl.pcolormesh(xygrid[:,0].reshape(shp), xygrid[:,1].reshape(shp), 
              np.exp(bg_ll).reshape(shp), cmap='Blues')

pl.xlim(0,1.)
pl.ylim(21,17.5)

The comparison!





In [ ]:



    


pl.figure(figsize=(4,8))

pl.pcolormesh(xygrid[:,0].reshape(shp), xygrid[:,1].reshape(shp), 
              np.exp(ll - bg_ll).reshape(shp), cmap='Blues')

pl.xlim(0,1.)
pl.ylim(21,17.5)

ND interpolation of isochrone





In [ ]:



    


_iso = iso[::-1]
tck, u = interpolate.splprep([_iso['gP1']-_iso['rP1'],_iso['gP1']-_iso['iP1'],_iso['gP1']-_iso['zP1']], 
                             u=_iso['iP1'], k=3, s=1E-4)
u_fine = np.linspace(u.min(), u.max(), 1024)
gr,gi,gz = interpolate.splev(u_fine, tck)

# tck, u = interpolate.splprep([_iso['gP1']-_iso['iP1']], 
#                              u=_iso['iP1'], k=3, s=1E-4)
# u_fine = np.linspace(u.min(), u.max(), 1024)
# gi = interpolate.splev(u_fine, tck)



    











In [ ]:



    


pl.figure(figsize=(4,8))

pl.scatter(itrp_gi, itrp_i, alpha=0.5, color='k')
pl.scatter(gz, u_fine+DM, alpha=0.5, color='r')
pl.scatter(gi, u_fine+DM, alpha=0.5, color='g')
pl.scatter(gr, u_fine+DM, alpha=0.5, color='b')

pl.xlim(0,1.)
pl.ylim(21,17.5)

Check multivariate gaussian





In [ ]:



    


bg_h = 0.05
bg_V = np.diag([bg_h]*2)**2
W = np.array([[1, -1],   # g-i
              [0, 1]])  # i magnitude

# each covariance C = WCW^T
bg_V = np.einsum('mj,jk->mk', W, bg_V)
bg_V = np.einsum('lk,mk->ml', W, bg_V)
bg_V



    











In [ ]:



    


_Cov = bg_Cov[:10]
_X = bg_X[:10]
_Cov.shape, _X.shape



    











In [ ]:



    


_X2 = xygrid
_Cov2 = np.zeros(_X2.shape + (_X2.shape[-1],))
_Cov2.shape, _X2.shape



    











In [ ]:



    


# METHOD 1:
_V = _Cov + bg_V[np.newaxis] # uncertainties plus smooth
_ll1 = log_multivariate_gaussian(_X[np.newaxis], _X2[:,np.newaxis], _V[np.newaxis])
_ll1 = logsumexp(_ll1, axis=1)
_ll1.shape



    











In [ ]:



    


# METHOD 2:
_V12 = _Cov[np.newaxis] + _Cov2[:,np.newaxis] + bg_V[np.newaxis,np.newaxis] # uncertainties plus smooth
_ll2 = log_multivariate_gaussian(_X[np.newaxis], _X2[:,np.newaxis], _V)
_ll2 = logsumexp(_ll2, axis=1)
_ll2.shape
# print(_V12.shape, _X.shape, _X2.shape, _ll2.shape)



    











In [ ]:



    


fig,axes = pl.subplots(1,2,figsize=(8,8),sharex=True,sharey=True)

axes[0].pcolormesh(xygrid[:,0].reshape(shp), xygrid[:,1].reshape(shp), 
                   np.exp(_ll1).reshape(shp), cmap='Blues')

axes[1].pcolormesh(xygrid[:,0].reshape(shp), xygrid[:,1].reshape(shp), 
                   np.exp(_ll2).reshape(shp), cmap='Blues')

pl.xlim(0,1.)
pl.ylim(21,17.5)



    











In [ ]:



    


def worker(allX, allCov, otherX, otherCov=None, smooth=None):
    if otherCov is not None:
        V = allCov[np.newaxis] + otherCov[:,np.newaxis]

        if smooth is not None:
            H = np.eye(allCov.shape[-1]) * smooth**2
            V += H[np.newaxis,np.newaxis]

        ll = log_multivariate_gaussian(allX[np.newaxis], otherX[:,np.newaxis], V)

    else:
        V = allCov
        if smooth is not None:
            H = np.eye(allCov.shape[-1]) * smooth**2
            V += H[np.newaxis]
        ll = log_multivariate_gaussian(allX[np.newaxis], otherX[:,np.newaxis], V[np.newaxis])

    ll = logsumexp(ll, axis=0) # NOTE: could also max here
    return ll



    











In [ ]:



    


worker(_X, _Cov, _X2).shape

Check isochrone func





In [ ]:



    


with h5py.File(XCov_filename, "r") as f:
    _X = f['isochrone']['X'][:]



    











In [ ]:



    


fig,axes = pl.subplots(1,4,figsize=(15,5), sharey=True)

DERPS = [_iso['gP1']-_iso['rP1'],_iso['gP1']-_iso['iP1'],_iso['gP1']-_iso['zP1'],_iso['rP1']-_iso['zP1']]

for i in range(1,4+1):
    axes[i-1].plot(_X[:,i], _X[:,0]+DM, ls='none', alpha=0.4)
    axes[i-1].plot(DERPS[i-1], _iso['iP1']+DM, ls='none', alpha=0.4)
    
axes[0].set_ylim(22,14)



    











In [ ]:

Content source: adrn/globber

Similar notebooks:

notebook.community | gallery | about